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ABSTRACT 

Here we present a simple, parameter-free, non-perturbative algorithm that gives low-redshift 
cosmological particle realizations accurate to few-Megaparsec scales, called MUSCLE (MUl- 
tiscale Spherical ColLapse Evolution). It has virtually the same cost as producing A^-body- 
simulation initial conditions, since it works with the ‘stretch’ parameter the Lagrangian 
divergence of the displacement field. It promises to be useful in quickly producing mock cat¬ 
alogs, and to simplify computationally intensive reconstructions of galaxy surveys. MUSCLE 
applies a spherical-collapse prescription on multiple Gaussian-smoothed scales. It achieves 
higher accuracy than perturbative schemes (Zel’dovich and 2LPT), and, by including the void- 
in-cloud process (voids in large-scale collapsing regions), solves problems with a single-scale 
spherical-collapse scheme. Slight further improvement is possible by mixing in the 2LPT es¬ 
timate on large scales. Additionally, we show the behavior of ^ for different morphologies 
(voids, walls, filaments, and haloes). A Python code to produce these realizations is available 
at http://skysrv.pha.jhu.edu/~neyrinck/muscle.html. 
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1 INTRODUCTION 

In the current paradigm, structures on extragalactic scales arise 
from the stretching and bunching-together of the ‘dark-matter 
sheet.’ A conceptual understanding of this process is a fundamen¬ 
tal goal of astronomy, and is closely tied to several other topics in 
cosmology and galaxy formation. A Lagrangian picture, following 
particles (i.e. mesh locations on this ‘dark-matter sheet’) (Shan- 
darin et al. 2012; Abel et al. 2012) instead of fixed comoving posi¬ 
tions, is essential to understanding this process. Approximations to 
full A-body dynamics are useful both for this understanding, and 
to produce approximate A-body realizations, e.g. for mock galaxy 
catalogs (e.g. Kitaura et al. 2014), and for fast explorations of large 
parameter spaces for Bayesian inference of the initial-conditions 
(IC) density field (e.g. Kitaura & EnBlin 2008; Kitaura et al. 2012; 
Kitaura 2013; HeB et al. 2013; Jasche & Wandelt 2013; Leclercq 
et al. 2015). 

The Zel’dovich approximation (Zel’dovich 1970, ZA), linear- 
order LPT (Lagrangian perturbation theory) already captures much 
of the essential physics of structure formation (e.g. White 2014), 
producing a cosmic web accurately to rather small scales. Improve¬ 
ments on ZA in various regimes have been proposed. Going to sec¬ 
ond order (2LPT) improves accuracy for small fluctuations, i.e. at 
large scales and early times, useful for IC generation (Scoccimarro 
1998; Crocce et al. 2006). At small scales, truncating the power 
spectrum shortward of the nonlinear scale (Kofman et al. 1992), 
or using the adhesion model (Kofman & Shandarin 1988; Kofman 
et al. 1992; Shandarin 2009; Valageas & Bernardeau 2011; Rid¬ 


ding et al. 2015) suppresses unphysical overcrossing of particles in 
collapsed structures. 

Many of these approaches work with the Lagrangian diver¬ 
gence of the displacement held, ^(q) = Vg • ^(q). Here, the dis¬ 
placement held is ^(q) = x(q) — q, where q denotes Lagrangian 
(IC) coordinates of a mass element, and x denotes its Eulerian, flnal 
coordinates. A recent approach (Neyrinck 2013, N13) works non- 
perturbatively with fj, using a low-density limit of the spherical- 
collapse (sc) evolution of a mass element, found by Bernardeau 
(1994). This SC relationship was also investigated by Protogeros & 
Scherrer (1997) and in the context of IC reconstructions by Mo- 
hay aee et al. (2006). 

Importantly, displacement-divergence (^-based) schemes giv¬ 
ing cosmological realizations are essentially as fast as producing 
the initial conditions for an A-body simulation. They have three 
steps: (1) Generate a pixelated linear-theory density held Sun at the 
desired redshift, consistent with a linear power spectrum. (2) Esti¬ 
mate ^ from S\in- (3) Generate the flnal displacement fleld ^ with 
an inverse-divergence operator, e.g. using an EFT; apply ^ to par¬ 
ticles on a regular lattice to get their flnal positions. 

In the ZA, = Sunl in 2LPT, there is a non-local 

functional giving roughly parabolic in the local 6\[n (N13). In 
the SC prescription, for each particle (IC density fleld pixel), 

V’sc(<5lin) = 3^1 - (2/3)<5lin - 3, <5lin < 3/2; 

= -3, <5,in ^ 3/2. (1) 

Or, (over)compactly, 

V’sc(<5iin) = 3 Rev/1 - (2/3)<5iin - 3. (2) 
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Setting ^sc = — 3 produces a collapsed volume element as closely 
as possible when dealing only with since ^ = —3 in three di¬ 
mensions if adjacent particles coincide. 

SC successfully reins in particles, preventing overcrossing on 
the interparticle scale, even when Sun is non-perturbatively large. It 
also gets densities remarkably right in voids, unlike ZA (overevacu¬ 
ating them), and 2LPT (underevacuating them) (Sahni & Shandarin 
1996). When applied far enough outside the perturbative regime 
of small fluctuations, 2LPT even produces overdensities in void 
centers (N13). However, when applied at the single, highest La- 
grangian resolution of a nonlinear-resolution density field, SC has 
problems as well. It gives realizations with low power on large 
scales; also, the Fourier-space cross-correlation coefficient departs 
from unity at larger scales than in ZA. 

Several other techniques have been proposed recently to pro¬ 
duce approximate, fast A^-body or halo density fields, such as 
PINOCCHIO (Monaco et al. 2002,2013), COLA (Tassev et al. 2013), 
and QPM (White et al. 2014). Kitaura & HeB (2013, KH13) found 
one insightful fix to the SC’s issues, calling it Augmented La- 
grangian Perturbation Theory (ALPT). In Fourier space, 

'0alpt(A:) = G(A:)'02lpt(/c) + [1 — G {k)]'ip sc (k), (3) 

an interpolation in scale between 2LPT, which excels on large 
scales with small fiuctuations, and SC, which gets small scales 
right. Here G{k) is a Gaussian smoothing kernel. ALPT is a key 
ingredient in the PATCHY algorithm (Kitaura et al. 2014), which 
produces mock galaxy catalogs using some other novel, useful pre¬ 
scriptions to treat redshift-space distortions and galaxy biasing. 


2 METHOD AND RESULTS 

Here, we give an alternative fix to the SC prescription: making it 
multiscale. As we show below, the trouble with the SC prescription 
is that it is applied at the single scale of the interparticle spacing, 
preventing voids in clouds (voids within larger-scale collapsing re¬ 
gions, e.g. Sheth & van de Weygaert 2004) from properly collaps¬ 
ing. 

MUltiscale Spherical-ColLapse Evolution (MUSCLE) addi¬ 
tionally checks scales larger than the interparticle scale for col¬ 
lapses, thus including the void-in-cloud process. Like the above 
schemes, MUSCLE works with step (2), the estimation of '0((5iin), 
with S\in defined on an interparticle separation (resolution of the 
initial density field), R[p. The prescription for -0 is 

^musc (^lin) — 3 1 ^ ^ Uud 

Gr((5iin) < 3/2, Vr ^ i^ip; 

= — 3, otherwise. (4) 

Gr((^iin) denotes the Sun field, Gaussian-filtered at scale r. Note 
that generally, 0musc will have nonzero mean, so an additional step 
is to subtract a constant from all V^musc > — 3 to give zero mean. 

In practice, a finite, and preferably small, number of r > R[p 
have to be tried. As implemented, we search for collapses at r = 

2'^Rip, for integers n ^ 0. The Gaussian smoothing is done in 
Fourier space, requiring an FFT and inverse FFT. We also tried 
using a cubic top-hat filter, averaging together sets of 8 particles 
on larger scales. This did visually localize halo regions better than 
the Gaussian smoothing, but the below cross-correlation performed 
more poorly. 


A large factor between smoothing scales searched for col¬ 
lapses gives greater speed, while a small factor captures more col¬ 
lapses, giving more accuracy. At our fiducial factor of 2, generating 
a random ACDM 256 h~^ Mpc realization with 256^ particles at 
z = 0 with the MUSCLE Python code searched 7 different smooth¬ 
ing scales; after no collapses were found at 64 h~^ Mpc, no larger 
scales were necessary. On one 2.6 GHz processor, generating this 
realization with ZA, 2LPT, and MUSCLE took 14, 28, and 59 sec¬ 
onds, respectively. 

Fig. 1 shows 2D patches of 3D particle realizations at redshift 
z = 0 using various approximate schemes from the initial condi¬ 
tions. It also shows the full A^-body results, and ‘Perfect 0’ - this is 
generated by measuring the actual 0 field from each particle’s N- 
body position using an FFT, and inserting it into step 2 above. This 
procedure zeroes the curl of the displacement field, but preserves 
all information in the displacement-divergence 0. See Chan (2014) 
for another investigation of zeroing the curl. In Fig. 2, a red line is 
drawn from each particle’s actual A^-body final position to that in 
the approximate scheme. 

In voids, densities (visible from particle separations) are too 
low in the ZA (Zeld), and too high in 2LPT. Haloes are puffy in 
Zeld, and even puffier in 2LPT. In this simulation, the interparticle 
scale is about 0.Sh~^ Mpc, with 256^ particles in a 200 h~^ Mpc 
box. The linear-theory variance on this scale much exceeds unity, 
^\in — 7.3; it is not surprising that a higher-order perturbative 
scheme (2LPT) fails worse than a lower-order scheme (Zeld) when 
the perturbative parameter is large. SC, as implemented in Eq. (1), 
gets densities right in voids, and pulls in overcrossed particle po¬ 
sitions in haloes, but as Eig. 2 shows, rather large-scale displace¬ 
ments are wrong in SC. In MUSCLE, the large-scale displacement 
problems are largely solved. This shows that the major deficiency 
in SC is in its treatment of the void-in-cloud process. 

We also show results of ‘MUSCLE-1-2LPT’, using the ALPT 
strategy of interpolating 0 between 02 lpt on large scales and 
0MUSCLE on small scales, as in Eq. 3. The best interpolation 
smoothing parameter using MUSCLE at = 0 seems to be about 
30 h~^ Mpc, judging by R(k) as shown below in Eig. 3. We also 
observed gains in R{k) in ALPT (using SC instead of MUSCLE) 
on small scales up to this smoothing length, although KH13 found 
that a smaller smoothing was optimal. So, we use the same smooth¬ 
ing length for each. Interpolating MUSCLE with 2LPT gives a few- 
percent boost in R(k) on non-linear scales, and gives perhaps a 
slight visual improvement over MUSCLE, in Pigs. 1 and 2. A smaller 
interpolation scale, <10/i“^Mpc, resulted in 2LPT visibly cor¬ 
rupting large haloes and voids. 

In Pig. 3, we show measures of the agreement between the 
A^-body result and various approximations. At top, we show the 
Pourier-space cross correlation, R(k) = PsxS' /VPsPs' between 
the z = 0 density field and the various approximate density fields, 
using nearest-gridpoint density assignment on a 128^ grid. PsxS' 
here is the cross-power between S and S'. Judging by R{k), sensi¬ 
tive to both Pourier amplitudes and phases, 2LPT performs a couple 
of per cent better than MUSCLE at /c ~ 0.2 h~^ Mpc. Thus, even 
though muscle’s detection of large-scale collapses somehow man¬ 
ages to treat the large-scale tidal field with some accuracy, 2LPT 
still seems to be the best available method for small fiuctuations, 
e.g. to generate high-redshift initial conditions. 

However, MUSCLE outperforms both perturbative schemes on 
small scales. Interpolating in scale between 2LPT and MUSCLE 
(‘MUSCLE-I-2LPT’) gives the highest R{k) at all scales. Interpo¬ 
lating between 2LPT and SC (labeled ‘ALPT’) gives results only 
slightly worse than interpolating with MUSCLE. Note that due to 
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Figure 1. Initially square 2D patches, 60/i“^Mpc on a side, of 
256^-particle realizations made using approximate schemes using the 
displacement-divergence In ‘MUSCLE+2LPT,’ V^muscle is interpo¬ 
lated in scale as in ALPT(shown here as well, using SC instead of MUS¬ 
CLE), both with smoothing lengths of 30 h~^ Mpc. ‘Perfect '0’ shows posi¬ 
tions after removing the curl of the displacement field, and shows ‘ A/^-body’ 
shows the actual simulation. 

cosmic variance, small differences in R{k) may not persist for all 
realizations (see KH13 for error estimates.) However, note that the 
R{k) curves may (co)vary together up and down more than inde¬ 
pendently, so small differences may be more statistically signifi¬ 
cant than the error bars suggest. The superiority of all Lagrangian 
schemes is obvious from the low ‘Eulerian linear theory’ curve 
(Tassev & Zaldarriaga 2012), essentially the propagator (Crocce 
& Scoccimarro 2006). 

The bottom panel of Fig. 3 shows the ratio of the nonlinear 



Figure 2. Same as Fig. 1, except with a red line joining the position of each 
particle in the full A/^-body result to the approximation. 


power spectrum to those in the various approximations. The power 
spectra in the t/;-based schemes differd only slightly, but the ap¬ 
proximation with power closest to full A^-body is ‘MUSCLE-I-2LPT.’ 
Notably, the large-scale power deficiency in SC is fixed in MUSCLE. 
In both measurements, while MUSCLE performs well, there is still 
much distance left to go to ‘Perfect 'ipj which would be the result 
if the final t/) could be perfectly predicted. 

MUSCLE also gets the 1-point PDF of the field far more 
accurately than the perturbative approaches. Fig. 4 shows mass- 
weighted histograms of densities computed using the Voronoi Tes¬ 
sellation Field Estimator (e.g. Schaap & van de Weygaert 2000), 
with each particle contributing once to the histogram. The agree¬ 
ment on the low-density tail between the A^-body realization and 
the SC and MUSCLE realizations is remarkable. Note even the un¬ 
physical low-density peak in 2LPT. At high densities, SC and MUS¬ 
CLE also perform the most accurately of the approximations, al¬ 
though they still fall well short of A/^-body. The MUSCLE PDF is 
the closest of the approximations, evidently since it captures addi¬ 
tional large-scale collapses compared to SC. But even ‘Perfect t/;’ 
falls somewhat short. Adding large-scale 2LPT modes to SC and 
MUSCLE gives almost no change in these PDF’s. 

Fig. 5 shows 2D histograms of both as predicted using 
MUSCLE, and as measured at = 0 in the simulation, versus its 
linear-theory prediction from the initial conditions, ^zeid- The re- 
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Figure 3. Top: Fourier-space cross-correlations of the A^-body density 
fields with the various approximations. Bottom: Ratios of power spectra of 
the various approximations to that of the full A/^-body density field. 
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Figure 5. A comparison of = Vq • ^ as predicted from the initial con¬ 
ditions using MUSCLE, to that measured in an A/^-body simulation. V^zeid^ 
on the x-axis of both plots, is the linear-theory (Zel’dovich) prediction. The 
colorbar numbers indicate the number of particles (out of 256^) in that two- 
dimensional bin. The white-and-black line indicates where the ZeTdovich 
prediction would lie. 


used elsewhere use the FFT method (which does not give a clear 
locus of points at t/; = —3, as seen here). 

Another interesting measurement is the distribution of t/; in the 
A/^-body simulation for different morphologies of particles: void, 
wall, filament and halo. To classify particles, we use the ORIGAMI 
algorithm (Falck et al. 2012), which assigns the morphology of a 
particle according to the number of orthogonal axes by which any 
other particle crosses it (void=0, wall=l, filament=2, and halo =3), 
comparing the initial and final conditions. 

Fig. 6 shows these PDF’s for the various morphologies. Except 
for ‘Halo,’ the distributions of ^ look quite Gaussian. One part of 
the ‘Halo’ distribution is a sharp spike at t/; = —3, characterizing 
a Lagrangian patch that precisely contracts to an Eulerian point. 
There is another component that looks somewhat Gaussian, as well; 
it would be interesting to differentiate these populations physically. 
One possibility is that the high-t/; tail consists of particles on their 
first infall. But also, the ‘Eilament’ population has a small bump at 
t/; = — 3, suggesting a small amount of contamination between the 
two morphologies as detected by ORIGAMI. 

The positions of the other curves are also informative. Erom 
void to halo, the mean t/;’s of each morphology are t/; = 

(1.80,0.58, —0.54, —2.25). If wall and filament volume elements 
typically collapsed along one and two axes, their non-collapsing 
dimensions staying fixed in comoving coordinates, their means 
would be at ^ — 1 and — —2. Their higher 'ip values indicate 
that their non-collapsing dimensions tend to stretch substantially. 


Figure 4. Mass-weighted density PDF’s of the various realizations, mea¬ 
sured using a Voronoi tessellation. The thick magenta dashdotted curve 
shows the PDF using raw SC, and the thick dashed and curve shows the 
PDF using MUSCLE. The nearly overlapping thin curves of the same colors 
show PDF’s in ALPT and ‘2LPT+MUSCLE,’ i.e. mixing in 2LPT on large 
scales. 


lationship is simple in MUSCLE: if pJzeid < —1.5, '0musc = —3. 
But there are also particles with 'ipzeid > —1.5 that get mapped 
to —3; for these, ^zeid dips below —1.5 when smoothing on some 
scale. There are two ways used in this paper to measure the diver¬ 
gence ip from the final conditions: with an FFT; and a ‘difference 
method,’ in which 'ip is measured at a particle by differencing parti¬ 
cle positions initially on either side of it (this has a resolution twice 
as coarse as the FFT method). The results shown here and in the 
next figure use the difference method, while the ‘Perfect t/;’ results 


3 DISCUSSION AND CONCLUSION 

Here we have presented a conceptually simple, non-perturbative 
Lagrangian scheme that produces fast A^-body realizations, essen¬ 
tially as fast as producing initial conditions of a simulation, and 
more accurately than other such schemes. It applies a spherical- 
collapse criterion to the Lagrangian divergence of the displacement 
field, 'ip, on the pixel scale of the initial conditions, and on larger 
scales as well. Approximate A^-body schemes have proven useful 
to produce mock galaxy catalogs (Kitaura et al. 2014). Some past 
and upcoming surveys include emission-line galaxies that do not 
necessarily occupy the largest dark-matter haloes, such as WiggleZ 
(Drinkwater et al. 2010) and eBOSS (Comparat et al. 2013). For 
these, fast approximations which accurately resolve the cosmic web 
to small scales will be increasingly important, for example to pro¬ 
duce mock galaxy catalogs with a biasing prescription (e.g. Kitaura 
et al. 2014; Ata et al. 2015). 

There is still room for improvement in such a t/;-based scheme; 
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Figure 6. Stretch parameters 'ip = Vg • ^ measured in an A/^-body sim¬ 
ulation of interparticle spacing ^ lh~^ Mpc, separated by morphological 
type. The distribution for each type looks rather Gaussian, except for ‘Halo;’ 
see discussion in text. 


if it were possible to perfectly predict ip, realization accuracies 
would increase substantially. One possible way forward is to pre¬ 
dict the final morphologies (void, wall, filament, halo) of particles, 
and assign 'ip to each differently, according to their values in Fig. 
6. Note that although SC and MUSCLE collapse haloes rather well, 
filaments do not visibly tighten substantially compared to the ZA, 
at the resolution presented here (at higher resolution, they would 
tighten some filaments, since MUSCLE limits overcrossing). How¬ 
ever, a separate prescription for each morphological type would 
sacrifice conceptual simplicity, and likely introduce some ad-hoc 
parameters. Another idea is to remap the approximate realization’s 
matter-density PDF to that of a full A^-body simulation (Leclercq 
et al. 2013); note that this step would be bypassed if generating 
a mock galaxy catalog. Also, we found that interpolating between 
2LPT and MUSCLE as in ALPT can slightly improve accuracy over 
MUSCLE, at the expense of slightly adding to the conceptual com¬ 
plexity and the run time. 

A Python package that generates all t/;-based particle realiza¬ 
tions discussed here, including an ALPT interpolation in scale, is 
available at http : / /skysrv. pha . jhu . edu/ ~neyrinck/ 
muscle . html. It interfaces with CAMB via the CosmoPy pack¬ 
age, at http : / /www. if a . hawaii . edu/cosmopy/. 
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